arXiv: 1503.04607vl [physics.comp-ph] 16 Mar 2015 


A Python Class for Higher-Dimensional 
Schrodinger Equations 

Amna Noreen, Member, lAENG, Kare Olaussen, Member, lAENG, 


Abstract —We announce a Python class for numerical solution 
of Schrodinger equations in one or more space dimensions, 
employing some recently developed general classes for numeri¬ 
cal solution of partial differential equations, and routines from 
numpy and scipy. sparse . linalg (or scipy. linalg for 
smaller problems). 

Index Terms —Quantum Mechanics, Anharmonic Oscillators, 
Sparse Scipy routines. 

I. Introduction 

URPRISINGLY many basic problems from physics and 
other natural sciences can be solved by one-dimensional 
analysis, mainly due to the symmetries of nature. A suc¬ 
cessful approach is often to search for symmetric solutions, 
or (through separation of variables for linear problems) 
solutions which can be composed of single-variable functions 
with simple symmetry behavior. 

However, not all problems of interest enjoy a high degree 
of symmetry. Many of them are still amenable to numerical 
analysis, although not always in any easily available manner. 
In this note we present an attempt to improve this situation 
for a particular class of problems from Quantum Mechanics, 
eigenvalue equations like 

[-A + V{r)]ipnir) = Eriipnir), (1) 

or variants and extensions of such equations. 

Our own interests in this class of problems arose when we 
attempted to generalize our method of very-high-precision 
solutions of such equations in one dimension Q, El, 0 
to higher dimensions. Our method works faster and most 
straightforward with some prior knowledge of eigenval¬ 
ues and eigenfunctions. For one-dimensional problems such 
knowledge can f.i. be obtained by use of the WKB method 
a, Q. The corresponding information is less accessible 
through analytic means in higher dimensions 0 , 0 . A 
numerical approach looks like a faster and more general 
strategy, which for most practical applications may anyway 
be sufficient. 

We believe that a numerical solver where the user can 
simply provide an essentially analytic formula for V{r) in 
Q, plus some information related boundary conditions and 
desired accuracy of the numerical approximation, would be 
of general interest. 

This can be realized using the freely available packages 
in numpy IS] and scipy ||9|, ifTOl . Three Python classes 
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developed to simplify the numerical solution of partial dif¬ 
ferentials in general have recently been announced M- 
We have refactored our original code to make use of these 
classes. 

For many problems of interesting complexity and size, the 
resulting code can be run on a normal laptop - and of course 
on more powerful computer systems. 

II. Basic numerical implementation 

Our code is constructed as a Python class, 
SchrodingerEquation which inherits the 
LatticeOperator class, and requires the Lattice 
and LatticeFunction classes, plus routines from the 
numpy and scipy packages. The latter packages are freely 
available for Windows, Mac OSX, and Linux platforms, 
or in (perhaps) numerically more efficient commercial 
versions. They allow a relatively simple formulation and 
organization of the problem in Python, with all heavy 
numerical calculations delegated to fast compiled routines. 

A. Lattice shape and geometry 

For discretization the continuum of positions r is replaced 
by a rectangular d-dimensional lattice of shape 

N =iNo,...,Nd-i), 

so that the total number of lattice points is TV = util Nk- 
It is our experience that a modern high-end laptop can 
handle models for which N < few x 10®. This allows 
resonably accurate treatment of one-particle problems in < 4 
dimensions, two-particle problems in < 2 dimensions, and 
three-particle problems in one dimension. 

Note that in Python all array positions are counted from 
zero. Hence each point of this lattice is described by an index 
vector 

n = {no,...,nd-i), where 0 < < iV^. 

Geometry can be introduced through a mapping n ^ Vn- 
The definition of a lattice model, along with a simple set 
of possible (rectangular) geometries, is implemented in the 
Lattice class. Two dim-dimensional parameters, rE and 
rO, can be given. They specify respectively the length of 
each lattice side (by default 1), and the position of the “lower 
left” corner of space (by default the origin). Note that all 
lattice points lie inside the spatial region defined in this way, 
half a lattice cell away from the edges. 

B. Boundary conditions 

Since the lattice must necessarily be finite for any finite- 
memory computer, it is necessary to define boundary condi¬ 
tions bC on -0 at its edges, or more precisely how tp should 



be extended beyond the edges. Some choices commonly used 
in Quantum Mechanics are: 

1) Periodic extension, 'P'. 

2) Symmetric extension, 'S'. 

3) Antisymmetric extension, 'A'. 

4) Extension with zero values, ' Z' . 

These choices may be different in different directions of the 
lattice, and — except for ' P' — different at the two edges 
of a given direction. All the above boundary conditions are 
implemented in the Lattice class. 

C. Lattice Laplacian. Stensil representations. 

The simplest implementation of a lattice Laplacian in one 
dimension is by use of the common discretization formula 

V [Xn) ~ -, (2) 

where A^; is the distance between nearest-neighbor lattice 
points. The generalizes to the standard {2d + l)-stensil 
definition of the lattice Laplacian, 

{ALlp){rr,) = 
d-1 

)] /drl (3) 

k=0 

Here dvk = dr [ k ] is the sidelength of the lattice cell in fc* 
direction, and 

^(-l-fc) = (^Ol ■ ■ ■ , tlk + 1, ■ ‘ ‘ ; tld-l), 

— (^0) ' ' ' 1 Xlk 1; ‘ ‘ ‘ 5 'ktd-l')- 

The discretization error of ([^ is of order 

d-1 

e^=0{Ydrl). (4) 

k=0 

With periodic boundary conditions equation Q is straight¬ 
forward to implement in numpy by use of the roll opera¬ 
tion. Other boundary conditions require more considerations 
and code, which have been incorporated in the Lattice 
and LatticeOperator classes. An arbitrary (short-range) 
operator O, represented by a stensil so{b) such that 

{Oip){rr,) =Y^o{b)'ipirn-b), ( 5 ) 

b 

is supported. Here ij^irn-b) is interpreted according to its 
boundary conditions when n — b falls outside the lattice. The 
stensil sO can be any d-dimensional numpy array provided 
by the user (but the program will run slowly if it is too large, 
since the sum over b is executed in native Python). By default 
sO is set to the stensil defined by equation ([^. 

III. Simple implementation 

It follows from the discussion above that a basic imple¬ 
mentation of the SchrodingerEquation class is quite 
straightforward. The lattice parameters can be set up by use 
of the Lattice class, and a standard approximation of the 
Laplace operator in equation Q is already available in the 
LatticeOperator class (or an alternative sensil can be 
provided). 


We further need to define an array representing the poten¬ 
tial L(r„) on every site of the lattice; this can be done in 
various ways by methods in the LatticeFunction class. 

We combine these to a method TplusVstensil (psi) 
in the class SchrodingerEquation. This performs the 
mapping 

'^{rn) ^-{ALil^){rn)+ V{rn)ip{rn) ( 6 ) 

for any given input array i/’(r„). This method can be used 
(i) to construct an explicit matrix representation of the 
linear operator in question, using the matrix () method, 
for further analysis by standard dense matrix routines in 
scipy. linalg, or (ii) as a LinearOperator alterna¬ 
tive to an explicit sparse matrix, for further analysis by the 
appropriate iterative routines in scipy. sparse . linalg. 
The implementation of the metod is straightforward and 
short: 

def TplusVstensil(self , psi) : 

"""Return (T + V) psi.""" 
psiO = -self .stensOp (psi) 
if self.def_V is not None: 

psiO += self . valuesV*psi 
return psiO 

The bulk of the code lies in the method stensOp, inher¬ 
ited from LatticeOperator, and methods in Lattice 
called by this routine. Note that all lattice dimensions and 
sizes (restricted by available computer memory), and all 
combinations of boundary conditions discussed above (ten 
possibilities for each direction), is handled transparently to 
the user. The stensil used to represent the operator T is 
stored in the dim-dimensional array self . stensil. It 
is by default the standard (2d -f l)-stensil for the lattice 
Laplacian, but is easily changed to any preferred array. 

A. Example: One-dimensional harmonic oscillator 

Consider the eigenvalue problem of the one-dimensional 
harmonic oscillator, 

- 'll:'.^{x) + ll)n{x) = Enlpnix). ( 7 ) 

The eigenvalues are £)„ = 2n -I- 1 for n = 0,1,..., and the 
extent of the wavefunction 4’n{x) can be estimated from the 
requirement that a classical particle of energy En is restricted 
to x"^ < En- A quantum particle requires a little more space. 

The numerical analysis of this problem is indicated by the 
code snippet 

myL = Lattice (shape= (2**7,), bC='Z', 
rE=(25, ), r0=(-12.5, )) 
defV = lambda x: x[0]**2 
myS = SchrodingerEquation (myL, 
def_V=defV) 

myS . varOp = myS . TplusVstensil 
H = myS. matrix () 
evals = eigvalsh(H) 

In lines 1 we define a geometric region of size 25, with 
origin in the middle, and cover it with a one-dimensional 
lattice of 2^ = 128 sites. In line 2 we define the potential 
V = x^, and use this and myL to define an instance of 
SchrodingerEquation in lines 3. The property varOp 
specifies a operator used by SchrodingerEquation in 
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many situations. This is set to TplusVstensil (different 
from the default) in line 4. A 2^ x 2^ dense matrix repre¬ 
sentation, H, of the operator is calculated in line 5. Finally 
all eigenvalues of H are computed in line 6. Here almost all 
computational work is executed in lines 5-6. 

The behavior of the resulting eigenvalues is shown in 
Fig.g 



Fig. 1. The 128 lowest eigenvalues of equation jTJ, computed with the 
standard 3-stensil approximation for the Laplace operator (here the kinetic 
energy T). The parameters are chosen to illustrate two typical effects: With 
bC=' Z' boundary conditions the hai'monic oscillator potential is effectively 
changed to V = oo for x > 12.5, thereby modifying the behavior 
of extended (highly exited) states. The effect of this is to increase the 
eigenenergies of such states, to a behavior more similar to a particle-in-box. 
This is visible for n > 80. The effect of using the 3-stensil approximation 
for T is to change the spectrum of this operator from to (the slower 
rising) (2/Aa;)^ sin^(fcAa;/2). This is visible in the sublinear rise of the 
spectrum for TV = 2^. 

For a better quantitative assessment we plot some energy 
differences, in Fig. 



Fig. 2. The discretization error of energy eigenvalues when using the 
standard 3-stensil approximation for the one-dimensional Laplace operator 
(here the kinetic energy E). There is no improvement in Fgo beyond 
certain TV, because the corresponding oscillator state is too large for the 
geometric region. For the other states the improvement is consistent with 
the expectation of an error proportional to A^. This predicts an accuracy 
improvement of magnitude 2^^ = 4096 when the number of lattice 
sites increases from TV = 2^ to TV = 2^^ for a fixed geometry. The 
eigenvalues are computed by the dense matrix routine eigvalsh from 
scipy. linalg. 

This brute force method leads to a dramatic increase in 
memory requirement with increasing lattice size. For a lattice 


with N — 2^ sites, the matrix requires storage of 4E double 
precision (8 byte) numbers. For m — 13 this corresponds 
to about ^ Gb of memory, for m = 14 about 2 Gb. The 
situation becomes even worse in higher dimensions. 

Assuming that we are only interested some of the lowest 
eigenvalues, an alternative approach is to calculate these by 
iterative routine eigsh from scipy. sparse . linalg. 
The most straightforward change is to replace lines 5-6 in 
the previous code with the snippet 

N = myL.size 

H = LinearOperator {(N,N), 

matvec^myS. linOp, dtype=float) 
evals = eigsh (H, which=' SM' ,tol=10*-^ (-8) , 
k=128 , return_eigenvect ors=False) 

which will compute the lowest 128 eigenvalues (still a rather 
generous amount). This allows extension to larger lattices, as 
shown in Fig. 



Fig. 3. The discretization error computed by the routine eigsh from 
scipy. sparse . linalg. For a fixed lattice size the discretization error 
is essentially the same as with dense matrix routines. However, with a 
memory requirement proportional to the lattice size (instead of its square) it 
becomes possible to go to much larger lattices. This figure also demonstrates 
(F 70 ) that the eiTor can be limited by boundary effects instead of a finite 
discretization length A^;. 

With a sparse eigenvalue solver the calculation becomes 
limited by available computation time, which often is a much 
weaker constraint. With proper planning and organization 
of calculations the relevant timescale is the time to analyze 
and publish results (i.e. weeks or months). The computation 
time is nevertheless of interest (it shouldn’t be years ifT^ l. 
We have measured the wall clock time used to perform the 
computations for Figs.|g[^ performed on a 2012 Mac Mini 
with 16 Gb of memory, and equipped with a parallellized 
scipy library. Hence, the eigvalsh and eigsh routines 
are running with 4 threads. The results are plotted in Fig. 

Here we have used the eigsh routine in the most straight¬ 
forward manner, using default settings for most parameters. 
This means, in particular, that the initial vector for the 
iteration (and the subsequent set of trial vectors) may not be 
chosen in a optimal manner for our category of problems. It 
is interesting to observe that eigsh works better for higher¬ 
dimensional problems. The (brief) scipy documentation 
O says that the underlying routines works best when 
computing eigenvalues of largest magnitude, which are of no 
physical interest for our type of problems. It is our experience 
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by the numerical approximation. For a numerical solution we 
only have to change lines 1-2 of the previous code snippet 
to 

1 myL = Lattice (shape= (2**7,)*2, bC='Z', 

rE=(25, )*2, r0=(-12.5, )*2) 

2 defV = lambda x: x[0]**2 + x[l]**2 

in two dimensions, and 

1 myL = Lattice (shape= (2**7,)*3, bC='Z', 

rE=(25, )*3, r0=(-12.5, )*3) 

2 defV = lambda x: x[0]**2 + x[l]**2 + 

x[2]**2 


Fig. 4. The wall clock time used to find the lowest 128 eigenvalues, for 
various systems and methods. We have also used the dense matrix routine 
eigvalsh to compute the eigenvalues of a 2^ X 2^ {N = 2^^) two- 
dimensional lattice; not unexpected it takes the same time as for a 2^“^ one¬ 
dimensional lattice. Somewhat surprisingly, with eigsh it is much faster to 
find the eigenvalues for two-dimensional lattice than for a one-dimensional 
with the same number of sites, and somewhat faster to find the eigenvalues 
for a three-dimensional lattice than for the two-dimensional with the same 
number of sites. 


that the suggested strategy, of using the shift-invert mode 
instead, does not work right out-of-the-box for problems of 
interesting size (i.e., where dense solvers cannot be used). 
We were surprised to observe that the computation time may 
decrease if the number of computed eigenvalues increases, 
cf. Fig.|5] 


Time to compute k eigenvalues 



Fig. 5. One may think that it takes longer to compute more eigenvalues. 
This is not always the case when the number of eigenvalues is small, as 
demonstrated by this figure. The default choice of eigsh is to compute 
fc = 6 eigenvalues. For our two- and three-dimensional problems this looks 
close to the optimal value, but it is too low for the one-dimensional problem. 


B. Example: 2- and 3-dimensional harmonic oscillators 
The (i-dimensional harmonic oscillator 

= En^nir), ( 8 ) 

has eigenvalues — {d 2n), for n = 0,1,.... The 
degeneracy of the energy level En is = (n + 1) in 
two dimensions, and = 4(n + l)(n + 2) in three 
dimension^] These degeneracies may be significantly broken 

*The general formula is . 


in three dimensions. 

As already discussed, the routine eigsh works somewhat 
faster in higher dimensions than in one dimension (for the 
same total number N of lattice points). The corresponding 
discretization error is shown in Figs. m 


Numerical errors in eigenvalues 



Fig. 6. The discretization error of energy eigenvalues when using the 
standard 5-stensil approximation for the two-dimensional Laplace operator. 
Exactly, the states and Ego are the two edges of a 13-member multiplet 
with energy 26, and the state £'12 is the middle member of a 5-member 
multiplet with energy 10. With the chosen parameters all states considered 
a well confined inside the geometric region; hence we do not observe any 
boundary correction effects. 


IQi 

10 ° 

10-4 

10-2 


Numerical errors in eigenvalues 


- I p(exact) p 
^ 1 

: A 

■ ► 





AAA 

E56 


► ►► 

Egs 


• •• 

Ei 5 

A 

TVT 

Eq 

► 



► 


log2 N ▼ 

J^I_l 


15 16 17 18 19 20 21 


Fig. 7. The discretization error of energy eigenvalues when using the 
standard 7-stensil approximation for the three-dimensional Laplace operator. 
Exactly, the states E^q and Es 3 are the two edges of a 28-member multiplet 
with energy 15, and the state E 15 is the middle member of a 10-member 
multiplet with energy 9. 

































The discretization error continues to scale like A^. This 
means that a reduction of this error by a factor 2^ = 4 
requires an increase in the number of lattice points by a 
factor 2'^ in d dimensions. This means that is becomes more 
urgent to use a better representation of the Laplace operator 
in higher dimensions. Fortunately, as we shall see in the next 
sections, better representations are available for our type of 
problems. 

IV. FFT CALCULATION OF THE LAPLACE OPERATOR 

One improvement is to use the reflection symmetry of each 
axis (x — )■ —X, y —)• —y, etc.) to reduce the size of the spatial 
domain. This reduces by a half, without changing the 
number of lattice points. 

A much more dramatic improvement is to use some 
variant of a fast fourier transform (fft); After a Fourier 
transformation, ■i/'(r) —>■ the Laplace operator turns 

into multiplication, (—Ai/;) (r) —)■ 'iA(fe)- This means that 
application of the Laplace operator can be represented by 
(i) a Fourier transform, followed by (ii) multiplication by 
fc^, and finally (iii) an inverse Fourier transform. Essentially 
the same procedure works for the related trigonometric 
transforms. 

These are also practical procedures for lattice approxima¬ 
tions, due to the existence of efficient and accurat^algoritms 
for discrete fourier and trigonometric transforms. The time to 
perform the above procedure is not significantly longer than 
the corresponding stensil operations. The benefit is that the 
Laplace operator becomes exact on the space of functions 
which can be represented by the modes included in the 
discrete transform. 

With the SchrodlngerEqpiation class it is easy to 
employ the FFT representation. This is built into the method 
TplusV, which is the default setting for varOp. Thus we 
only have to comment out (or remove) line 4 in the first code 
snippet above: 

4 # myS.varOp = myS.TplusVstensil 

The obtainable accuracy with this procedure increases 
dramatically, as illustrated in Figs. [8pT| 

It might be that the harmonic oscillator systems are partic- 
ulary favorable for application of the FFT representation. One 
important feature is that the fourier components of the har¬ 
monic oscillator wavefunctions vanishes exponentially fast, 
like e“^ with increasing wavenumbers k^. This feature 
is shared with all eigenfunctions of polynomial potential 
Schrodinger equations, but usually with different powers of 
fc in the exponent (which may lead to a quantitative different 
behavior). 

Further, for systems with singular wavefunctions the corre¬ 
sponding fourier components may vanish only algebraically 
with fc^. The dramatic increase in accuracy cannot be ex¬ 
pected for such cases. 

^The error of a back-and-forth FET is a few times the numerical accuracy, 
i.e. in the range 10“^'* — 10“^® with double precision numbers. However, 
when an error of this order is multiplied by it can be amplified by several 
orders of magnitude. Hence, the range of fc^-values should not be chosen 
significantly larger than required to represent t/>(r) to sufficient accuracy. 


Numerical errors in eigenvalues 



Fig. 8. With a FFT representation of the Laplace operator the discretization 
error drops exceptionally fast with Aa; oc When it becomes “small 

enough” the effect of numerical roundoff becomes visible; the latter leads 
to an increase in error with A^. The results in this figure is for a one¬ 
dimensional lattice, but the behavior is the same in all dimensions. The 
lesson is that we should make Aa: “small enough” (which in general may 
be difficult to determine a priori), but not smaller. It may be possible to 
rewrite the eigenvalue problem to a form with less amplification of roundoff 
errors. 


Numerical errors in eigenvalues 



Eig. 9. Accuracy of computed eigenvalues for a ID oscillator. This figure 
may suggest that an increase in the number of lattice size N will lead to a 
accurate treatment of states with higher n. Our findings are that this is not 
the case: The results for N = 2^ and iV = 2® have essentially the same 
behavior as for V = 2®. 

V. Anharmonic oscillators 

The SchrodingerEquation class works with any 
computable potential. All one has to do is to change the 
definition of the function assigned to def_V. However, in 
most cases we no longer know the exact answer; this makes 
it difficult to assess the accuracy of the result. 

One simple test is to consider two-dimensional potentials 
of the form 

V{x,y)=x‘^+ cx^y'^ ^y"^, (9) 

for various values of c. For c = 0 and c = 6 this system 
an be separated to a set of two one-dimensional anharmonic 
oscillators (with x'^ type potential), and for c = 2 it can 
be separated in cylinder coordinates. Here comparison with 
essentially exact solutions of the separated one-dimensional 
problems can be made. Even in the absence of these, one 
may check for degeneracies in the spectrum. 

















Numerical errors in eigenvalues 



Fig. 10. Accuracy of computed eigenvalues for a 2D oscillator. As can 
be seen, a large number of the lowest eigenvalues can be computed to an 
absolute accuracy in the range with lattice of size 2® x 2®. 

We observe not improvement by going to 2^ x 2^ lattice, but no harm either 
(except for an increase in the wall clock execution time from about 3 to 30 
seconds for each combination of boundary conditions). 


Numerical errors in eigenvalues 



Fig. 11. Accuracy of computed eigenvalues for a 3D oscillator. As can 
be seen, a large number of the lowest eigenvalues can be computed to an 
absolute accuracy in the range 10~^^-10~^^ with lattice of size 2® x 2® x 
2®. We observe no improvement by going to 2^ x 2^ x 2^ lattice, but no 
harm either (except for an increase in the wall clock execution time from 
about 6 to 95 minutes for each combination of boundary conditions). 


VI. Boundary conditions for radial operators 


For some problems one may perform a partial (or full) 
symmetry reduction of equation 0. An example is 

(dp Id \ f‘ 

“(x:2 + zx: + t:2)+Z2 + 


r dr dz^ 


'4’{r,z) 

= E^{r,z). (10) 


In this case there is a natural boundary at r = 0. This is also 
a singular line for the equation. 

What is the natural boundary condition at r = 0? Equation 
(Hoj is often symmetric under r —> —r, depending on 
the form of V. With this symmetry equation ( [TOl i may be 
extended to r < 0, and the solutions classified according to 
their transformation under r —^ —r. This makes symmetric 
or antisymmetric boundary conditions the namral choice. 

Surprisingly, discussions of numerical approximation 
schemes for (very common) singular equations like ( [TOl i are 
difficult to find in the textbook literature. It may be tempting 


to introduce a new wavefunction, ij;{r,z) = ip{r, z), 

in order to eliminate the first order derivative in ( [Tol l. This 
would transform the equation into the general form 0. 
However, this would also make the searched-for solution 
(p(r^ z) singular at r = 0, hence difficult to approximate 
numerically. 

VII. Domains of general shape 

There is also a convenient way to define domains of 
general shape in numpy, by specifying a boolan vector which 
is True for all points in the domain, and False for all 
points outside. In this case we find only the ' Z' boundary 
conditions to be a general and unambiguous option. Stensil 
approximations of the Laplace operator may be the best 
choice for such cases. 
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